Numerically exact Green functions from Hirsch-Fye quantum Monte Carlo simulations 
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We present a new method for extracting numerically exact imaginary-time Green functions from 
standard Hirsch-Fye quantum Monte Carlo (HF-QMC) simulations within dynamical mean-field 
theory. By analytic continuation, angular resolved spectra are obtained without the discretization 
bias previously associated with HF-QMC results. The method is shown to be accurate even at very 
low temperatures (T = W/8QQ for bandwidth W) in the strongly correlated regime. 
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Photoemission spectroscopy (PES) and related tech- 
niques (inverse PES, X-ray absorption spectroscopy) are 
among the most useful experimental probes of elec- 
tronic properties of solids Already the angle- and 
spin-integrated variants can characterize materials as 
weakly or strongly correlated metals or as insulators. 
In addition, angle-resolved photoemission spectroscopy 
(ARPES) gives a nearly direct view on the occupied part 
of the electronic band structure. However, the inter- 
pretation of such experimental data and the separation 
of surface from bulk contributions is by no means triv- 
ial. Thus, both a reliable analysis of the experiments 
and an understanding of the underlying physics require 
comparisons with theoretical predictions. Often, good 
agreement is found with spectra obtained within den- 
sity functional theory (DFT), e.g., within the local den- 
sity approximation (LDA) [4]. For strongly correlated 
materials, however, the LDA and/or the interpretation 
of the DFT parameters as many-body dispersion break 
down; then, true many-body approaches, usually based 
on Hubbard-type models, are required. 

A nonperturbative treatment of Hubbard-type mod- 
els for correlated electron systems is possible by iter- 
ative solution of the DMFT self-consistency equations 
[§| using quantum Monte Carlo (QMC) methods. In 
the Hirsch-Fye QMC algorithm [J], the imaginary-time 
path integral is discretized into A time slices of uniform 
width At = /3/A; a Hubbard- Stratonovich (HS) transfor- 
mation replaces the electron-electron interaction at each 
time step by a binary auxiliary field which is sampled by 
standard Markov Monte Carlo techniques. After conver- 
gency, estimates of the local imaginary-time Green func- 
tion G{t) are obtained (only) on the grid t; = /Ar with 
(0 < Z < A). For the computation of spectral functions, 
this data has to be continued to the real axis, typically 
using maximum entropy methods (MEM) [5]; these reg- 
ularize the ill-conditioned inversion of the equation 
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by finding a spectrum A{u)) that is as smooth (or sim- 
ilar to a chosen default model) as possible within the 
constraints set by the data {G(t;)} and its error bars. 



If needed, the real part of the local Green function 
is then obtained by Kramers-Kronig transformation of 
ImG'(ti;) = —■kA{u!). Finally, momentum resolved spectra 
corresponding to ARPES measurements may be obtained 
via the real-frequency self-energy S(a;). 

While impressive results have been obtained using the 
QMC/MEM procedure as outhned above, e.g., in the 
context of "kink" anomalies in ARPES spectra it 
has one fundamental problem: the Hirsch-Fye QMC esti- 
mates of the imaginary-time Green function contain sys- 
tematic Trotter errors which are typically much larger 
than the statistical errors, often by orders of magnitude. 
However, the MEM takes only the statistical errors into 
account, partially by quite elaborate formalism, while the 
larger systematic errors in the data are neglected. Con- 
sequently, all resulting spectra are biased, to an unknown 
extent, by the Trotter discretization. 

In this Letter, a method is proposed which essen- 
tially eliminates these problems: using a novel extrap- 
olation scheme, we extract continuous estimates of the 
imaginary-time Green function without significant Trot- 
ter error from conventional HF-QMC data; these are then 
proper bases for analytic continuation techniques. The 
method works well even at very low temperatures and 
for comparatively coarse imaginary-time discretizations. 
These properties make it very attractive for future calcu- 
lations of spectra, e.g., in ah initio LDA-I-DMFT studies; 
they also establish that - contrary to common belief - 
reliable HF-QMC results do not depend on a good reso- 
lution of the rapid initial decay of G(t). 

Raw HF-QMC results - In the following, the method 
will be defined and illustrated using a quite ambitious ex- 
ample: the half-filled single-band Hubbard model (semi- 
elliptic density of states with bandwidth VF = 4) at the 
very low temperature T — 1/200 for U = 4.95, i.e., in the 
strongly correlated metallic regime. As HF-QMC results 
had so far only been reported for temperatures T > 1/50, 
these parameters have been believed to be out of reach 
of the Hirsch-Fye QMC method 0]. 

Discrete HF-QMC estimates of the imaginary-time 
Green function are shown as symbols in the main panel 
of Fig. [1] for relatively large values of the discretization 
At G [0.4,1]. Clearly, the Trotter errors in these data 
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FIG. 1: (Color online) Imaginary-time Green function for 
T = 1/200, U = 4.95: HF-QMC estimates for different 
discretizations Ar (symbols); lines are guides to the eye 
only. Also shown; Green functions corresponding to a model 
self-energy (solid line) and to the noninteracting limit (thin 
dashed line); inset: same data up to symmetry point r = (3/2. 

sets are very significant at least for t < 5, much larger 
than the statistical and convergency errors (of about 
10^^). However, extrapolation schemes that have suc- 
cessfully been used for observables such as quasiparticle 
weight Z, energies, and double occupancies [1, W\ can- 
not directly be applied here, since the r grid is different 
for each data set. In addition, none of the data sets re- 
ally covers the initial rapid drop of \G{t)\: even for the 
finest discretization Ar = 0.4, |G'(r)| is already halved 
at the first nontrivial data point. In fact, according to 
the assumption Q that a good resolution of the highly 
curved initial region was necessary, HF-QMC results for 
such coarse grids should not yield useful information at 
all; instead one would have to resort to much finer grids 
of At<1/(5[/)«0.04 [iI. 

These considerations, however, overlook an important 
physical point: The behavior of G{t) at small r is re- 
stricted by second-order weak-coupling perturbation the- 
ory (which is exact for the curvature of G at r ^ in 
the half-filled case; see below and Fig. [3]). Consequently, 
reasonable accuracy at small r can be expected from 
a model Green function based on a (here particle-hole 
symmetric) two-pole approximation to the self-energy, 
5^modei(w) = —Au)/ (w^ - Wg); since the weight A = U'^ 1^ 
is fixed by the leading large-frequency asymptotics [lOj, 
only the position of the poles (at iwo) remains as a free 
parameter. As indicated by the solid line in Fig. [TJ this 
simple model (for = 1-3) has all the features that we 
expect from the true Green function, including the rapid 
initial drop; in fact, a visual inspection of the trends sug- 
gests that it better approximates the true Green function 
(for small r) than any of the QMC data sets. Still, as we 
will show in the following, the latter contain the neces- 
sary information for computing very precise estimates of 
the true Green function at all r G [0,/?]; the model will 
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FIG. 2: (Color online) Green function G(r) for T = 1/200, 
U = 4.95: HF-QMC data of Fig. [1] (symbols) with interpo- 
lating fits (dashed/dotted lines) and numerically exact result 
after extrapolation Ar (thick solid line). Inset: Green 
functions relative to Gat, thin solid line: CT-QMC result 0]. 

only be needed for interpolating the raw QMC data. 

Extrapolation procedure - As a starting point for the 
extrapolation, we need accurate QMC data {G'(r/)} with 
reliable error bars. While arithmetic averages are clearly 
appropriate for combining the Green functions obtained 
in parallel QMC runs for the same impurity model, errors 
can be minimized (most notably in insulating phases) 
by using geometric averages for combining estimates of 
{G(r/)} from different DMFT iterations. In practice, we 
transform all data to log[— G(r)] before taking arithmetic 
averages; the appropriateness of this logarithmic scale 
will become apparent below (cf. Fig. [4]). 

In a second step, the averaged QMC data sets (for dif- 
ferent Ar) need to be transformed to a common grid. 
For this purpose, optimized two-pole approximations (as 
specified above) are obtained for each data set. Each dif- 
ference Gqmc ~ Ginodci (defined only on the specific grid 
{ZArj^p) is interpolated by a natural cubic spline which 
is evaluated on a fixed, much finer grid; finally, Gmodei(T) 
is added to these intermediate results . As seen in Fig. 
[21 these interpolated curves (dashed and dotted lines in 
the main panel) are smooth and vary systematically, with 
a discretization dependence which vanishes at small r. 

In the third step and final step, least squares fits are 
performed for log[— G(r)] at each value of r (on the fine 
grid), taking quadratic and quartic contributions in Ar 
into account 'l3'|. The result of this procedure is shown 
as thick solid line in Fig. [21 In the inset, the approximate 
low-r asymptotics of the Green function, 

Gfit,(r) = -0.5 cxp [-2.36 r (1 + 0.36 r + 0.16 r^)-!] 

have been subtracted. At this scale, the discretization 
error in the QMC data is clearly visible; in contrast, our 
extrapolated result (thick solid line) is hardly distinguish- 
able from a continuous-time QMC estimate 7] (thin solid 
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FIG. 3: (Color online) Second and first order derivatives 
(main panel/ inset) of the Green functions of Fig. (2] 
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FIG. 4; (Color online) HF-QMC estimates of imaginary-time 
Green functions for T = 1/45, U = 5 for tiie metallic (upper 
set of curves) and insulating (lower set) phases; side bands 
(thinner lines) indicate corresponding error estimates. 



line). The fluctuations (on a much larger r scale for HF- 
QMC than for CT-QMC) suggest that both results have 
similar precision. The competitiveness of HF-QMC (plus 
extrapolation) at this extremely low temperature, i.e., for 
rather coarse At grids, may appear surprising. However, 
as shown in the inset of Fig. [21 the strongest absolute de- 
viations of the raw HF-QMC data occur at r ^ 1 (while 
the relative errors are largest at r « 2), i.e., at a position 
independent of At. This establishes that the rapid ini- 
tial decay of G(r) does not set the scale for useful values 
of At in HF-QMC and explains the good performance of 
our extrapolation procedure. 

The uniform convergence of the (interpolated) HF- 
QMC Green functions, in turn, may be traced back to 
the fact that their curvatures are asymptotically exact 
for T 0. Indeed, the second derivatives cP'G{t) / cLt'^ are 
visibly At dependent in the main panel of Fig. [3] only for 
T « 0.5; they all approach the correct asymptotic limit 
(PG{T)/dT'^\r=o+ = (1 + C/2/4)/2 and agree at large t 
within statistical errors. In contrast, as shown in the in- 
set of Fig. [3l the Trotter errors of the first derivatives 
do not vanish at t ^ 0; concurrently, the asymptotic 
exact value dG{T) I dT\T=o-\. is nonuniversal, i.e., depen- 
dent on temperature and phase (metallic or insulating). 
Note that the extrapolated HF-QMC results in Fig. [3] 
(thick solid lines) are hardly distinguishable from the cor- 
responding CT-QMC results (thin solid lines). 

The specific form of the Green function extrapolation 
procedure detailed above is based on the insight that the 
Green function varies on a logarithmic scale. This is 
clearly seen in Fig. |4] for a moderately low temperature 
T = 1/45 and similarly strong interaction U = h (close 
to the thermodynamic Mott transition) : While all curves 
become indistinguishable for t — > 0, both differences be- 
tween metal and insulator and the At dependence in the 
insulating phase involve many orders of magnitude; in 
the latter case, even the error bars span nearly an order 
of magnitude. Obviously, a direct extrapolation At — > 
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FIG. 5: (Color online) Green functions: absolute deviations 
of finite-Ar data of Fig. [4] from the extrapolated results for 
insulating/metallic phase (upper/lower panel). 



of the insulating Green function (instead of log[— G(t)]) 
would be hopeless for t > 5. Note that our choice guaran- 
tees the correct sign for the extrapolated Green function. 

As seen in Fig. [51 the absolute Trotter errors (i.e., the 
deviations at finite At from the numerically exact re- 
sults) of all data sets peak, again, at the same position 
(here t w 1.2, indicated by vertical lines). The remark- 
able constancy of this peak position, even across phase 
transitions, implies that the extrapolation technique can 
be trusted to capture the 1ow-t features of the Green 
function without bias, even for relatively coarse grids At. 
In comparison, the region t w /3/2 is more challenging (in 
the insulating phase): since HF-QMC results for t > 0.4 
(dotted line in Fig. [4]) are too metallic, a good resolution 
of the exponential decay of G(t) (see lower solid line in 
Fig. [4]) can only be expected from extrapolations which 
include small discretizations At < 0.3. 

Spectra on the real axis - So far, we have only specified 
how to extract Green functions in the imaginary-time do- 
main. The analytic continuation of this data to the real 
axis is still a highly nontrivial task; in fact, the extrap- 
olated imaginary-time data arguably deserves even more 
careful and sophisticated continuation procedures than 
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FIG. 6: (Color online) Main panel: Local spectral func- 
tion ^(a;) for T = 1/200, U = 4.95 (solid line) in compar- 
ison with noninteracting spectrum (dashed line). Left inset: 
momentum-resolved spectra j4(k,a;) in same energy range for 
momenta k = F to k = X. Right inset: same data as in inset 
of Fig. [5] except for the thin solid line which corresponds here 
to the spectrum of the main panel via Eq. [1] 



regular HF-QMC results, since - for the first time - the 
Green functions are exact within error bars. Adapted 
maximum entropy procedures and methods for efficient 
error analysis will be discussed elsewhere . In the case 
of very precise data (as in our examples), however, one 
may also neglect the small statistical errors and use Fade 
methods. Specifically, the local spectral function shown 
in the main panel of Fig. [S] has been obtained via a Fade 
approximant 15] to the imaginary- frequency self-energy 
Yi{iujn) using the first 199 positive Matsubara frequen- 
cies. This procedure is stable [as has been checked via 
sum rules for S(a;) and G{lij) and by varying the Fade 
parameters] and asymptotically exact in the noninter- 
acting limit. In addition, it gives direct access also to 
momentum-resolved spectra (corresponding to ARFES 
measurements), as visualized in the left inset of Fig.[6l 

The expected pinning of the quasiparticle peak height 
in A{lij) to its noninteracting value (cf. dashed line in 
main panel) is accurately observed; the spectrum is also 
reasonably smooth without any indications of artefacts. 
A(k, bj) clearly resolves the dispersions of the heavy 
quasiparticles (with m* /ni « 9) and of the incoherent 
Hubbard bands. So the results appear reasonable - but 
are they significantly better that those attainable using 
conventional HF-QMC /MEM procedures? This is indeed 
the case, as demonstrated in the right hand inset of Fig. 
[6l here the thin solid line, computed from the final spec- 
trum via Eq. [U is hardly distinguishable from the nu- 
merically exact G{t) (thick line; cf. Fig. ^ while the 
deviations of finite-Ar data (on which spectra would be 
based in conventional methods) are larger by several or- 
ders of magnitude. Thus, the spectrum of Fig. [6] is seen 
to be unbiased and numerically exact ~ for a temperature 
which had been deemed out of reach of HF-QMC. 

Discussion - We have presented a method that, fi- 



nally, allows to obtain numerically exact Green func- 
tions and (k resolved) spectral functions from regular 
Hirsch-Fye QMC calculations. Due to the simplifica- 
tions in DMFT, transport properties such as the opti- 
cal conductivity a{uj) could be computed without ad- 
ditional effort. Since the method has been successfully 
tested also for doped systems [IJ], it should greatly im- 
prove the attainable quality of spectra, in particular in 
the context of LDA+DMFT calculations. In principle, 
the recently developed continuous-time QMC methods 
should yield spectra of similar quality with comparable 
effort; however, this might require regularization pro- 
cedures which fail in insulating phases [l6| . Numeri- 
cal renormalization group methods appear hampered by 
their coarse high-frequency resolution and systematic er- 
rors. It contrast, detailed comparisons of numerically ex- 
act HF-QMC spectra with ground state estimates from 
dynamical density-matrix renormalization group are ex- 
pected to lead to new physical insight. Specificly, one 
may hope to thereby resolve the still controversial ques- 
tion of how the structure of the Hubbard bands changes 
across the Mott transition. 

Stimulating discussions with P.G.J, van Dongen and 
support by the DFG within the Collaborative Research 
Centre SFB/TR 49 are gratefully acknowledged. 
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